
function dxx = Cr3bp_Eqm2D(tau , xx , mu)
%
% CR3BP动力学模型，仅状态，2D
%
% 输入：
% t:               [1x1]            时间(s)
% xx:             [6x1]            状态(km , km/s)
% mu:           [1x1]            引力常数 (km^3/s^2)
%
% 输出：
% dxxdt:        [6x1]            动力学方程(km , km/s^2)
%
% 例：
% xx0 = [0.75 , 0 , 0.15 , 0];
% options = odeset('RelTol',1e-8,'AbsTol', 1e-8);
% [tt , xx] = ode113(@crtbpEqm2D , [0 , 10] , xx0 , options , aux);
% 
% 时间：2024年06月03日
% 作者：张晨
% 单位：CSU
% 邮箱：chenzhang@csu.ac.cn
%%%%%%%%%%%%%%%%%%%%%%%

% % 积分时间太长，跳过
% if nargin == 4
%     ode_t = clock;
%     if etime(ode_t , ode_t0) > 1
%         dxx = [];
%         return
%     end
% end

% 提取变量
x = xx(1);
y = xx(2);
vx = xx(3);
vy = xx(4);

% % 配置参数
% mu = aux.ems.mu;
% LU = aux.ems.LU;
% P1_req = aux.planet.req(3);
% P2_req = aux.planet.req(10);

% % P1/P2碰撞检测
% if aux.impact_IO == 1
%     impactP1 = sqrt( (x + mu)^2 + y^2 ) - P1_req / LU;
%     impactP2 = sqrt( (x + mu - 1)^2 + y^2 ) - P2_req / LU;
%     if impactP1 * impactP2 < 0
%         dxx = [];
%         return;
%     end
% end

% 状态动力学方程
r1_3pow = ((x + mu)^2 + y^2)^(1.5);
r2_3pow = ((x + mu - 1)^2 + y^2)^(1.5);
dxx = [vx;
    vy;
    x + 2 * vy - (1 - mu) * (x + mu) / r1_3pow - mu * (x + mu - 1) / r2_3pow;
    y - 2 * vx - (1 - mu) * y / r1_3pow - mu * y / r2_3pow];

end
